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I present an overview of pulse propagation methods used in nonlinear optics, covering both full- 
field and envelope-and-carrier methods. Both wideband and narrowband cases are discussed. Three 
basic forms are considered - those based on (a) Maxwell's equations, (b) directional fields, and (c) the 
second order wave equation. While Maxwell's equations simulators are the most general, directional 
field methods can give significant computational and conceptual advantages. Factorizations of the 
second order wave equation complete the set by being the simplest to understand. One important 
conclusion is that that envelope methods based on forward-only directional field propagation has 
made the traditional envelope methods (such as the SVEA, and extensions) based on the second 
order wave equation utterly redundant. 



Contents 



VIII. Acknowledgments 



I. Introduction 

II. Fields vs Envelopes 

A. The definition 

B. The big advantage 

C. Nonlinear polarization terms 

D. Multiple carriers, wideband envelopes 

E. Directionality 

F. Moving frames 

G. Estimating the computational cost 

H. Disadvantages 

III. Maxwell's equations 

A. Envelopes 

B. Transverse effects 

IV. Directional Maxwell's equations 



1. Special case: x 

A. Envelopes 

B. Transverse effects 



(2) 



V. 



Second order wave equations 

A. Traditional approach 

B. Time propagated direct solution 

C. Short pulse equation (SPE) 

D. Factorization approach 

1. Simple factorization 

2. Improved factorization 

3. Linear factorization 

4. Special case: x^^^ 

5. Factorization and envelopes 

6. Factorized fields 

E. Transverse effects 

VI. Forward-backward coupling 
VII. Conclusions 



3 
3 
3 
1 



i 



i 
i 
i 
i 

i 

ill 
ill 
iil 

1 

B 



'Electronic address: ' Dr . Paul. Kinsler@phy sics . org| 



References 



I. INTRODUCTION 

Here I discuss three important ways to tackle pulse 
propagation in nonlinear optics. These include methods 
that both do and do not follow the traditional approach 
of using pulse envelopes. The description is taken in the 
ID limit, but some discussion on including transverse 
effects is made. The aim is to cover the considerations 
relevant when modeling wideband fields, a regime not 
comprehensively treated in many standard texts 0, [3, 

The three ways are solving (a) Maxwell's equations, 
(b) directional Maxwell's equations, or the (c) standard 
second order wave equation. Solving Maxwell's equa- 
tions is a well established approach, with a long his- 
tory (i.e. finite difference time domain (FDTD), see 
6-g- |21)j although it is computationally intensive and 
has generally been little used in nonlinear optics (but 
see e.g. [1, ll^l)- Practical versions of directional 
Maxwell's equations have appeared only recently, such 
as that of Kolesik et al. [Til, [13] ; other approaches fol- 
lowed [i3, [i3) [ill- However the first proposal dates 
back to Fleck in 1970 [i^, although only as something 
of a remark in passing, rather than a full investigation. 
The most common approaches nonlinear optics are those 
based on the standard second order wave equation, par- 
ticularly with regard to envelope propagation and the cel- 
ebrated slowly varying envelope approximation (SVEA). 
The SVEA allows us to convert the second order wave 
equation into a first order equation that can efficiently 
propagate narrowband pulses. Recently the SVEA has 
been relaxed fl3, [il, Il9l |. extending the use to moderate 
bandwidths. However, much better approaches based on 
factorizing the second order wave equation also exist. An 
early example can be seen in but also most notably 
by Blow and Wood 2C|, and also the recent Ferrando et 
al. [lil and Genty et al. [H]. 

We can try to solve any of these equations directly. 
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without recourse to an envelope and carrier representa- 
tion. This means ensuring sufficient numerical resolution 
to integrate each of the field oscillations as it propagates 
across the simulation window. This approach is the stan- 
dard one when solving Maxwell's equations (i.e. FDTD), 
but generally in nonlinear optics an envelope approach 
is used. This has a number of advantages: it imposes a 
direction on the modeled pulse, and it removes the fast 
oscillations at the centre frequency. In combination with 
a moving frame, it can turn a pulse of rapidly oscillat- 
ing fields moving at the speed of light into a smooth, 
nearly-stationary waveform - with commensurate gains 
in simulation speed. These benefits usually come with a 
restriction on the allowed bandwidth of the pulse being 
modeled. 

This paper is organized as follows: in section [ill I com- 
pare field and envelope approaches. In section [TTTl I con- 
sider Maxwell's equations in both field and envelope pic- 
tures, followed in section IIVI by the same, but utilizing 
a directional rewriting of Maxwell's equations. Then, in 
section |V] I consider the role of second order wave equa- 
tions, in particular using factorization methods. Finally, 
in section IVTIl I present some conclusions. 

Although not directly relevant to the discussion here, 
it is also worth noting that directional waves in unstable 
resonators were quantized by Brown and Dalton(23j. 



II. FIELDS VS ENVELOPES 

It is often remarked upon that envelope methods work 
surprisingly well. However, this surprise seems to be 
based largely upon the SVEA equation for pulse propa- 
gation, which indeed contains many approximations (see 
18, 24 1). Recently it has been shown by several groups 
11,1s, [1^, [2l| that equations nearly identical to those 



generated by the SVEA can be found by assuming lit- 
tle more than the lack of backward-going field compo- 
nents. Even when revisiting Blow and Wood [1^, we 
can see that their mathematics contained minimal con- 
straints on the bandwidth of the envelope, although their 
specific nonlinearity model did contain such restrictions. 



A. The definition 

For envelope methods, the direction is imposed by the 
form for the carrier function, and is usually a plane wave 
traveling in the chosen direction. Thus the typical enve- 
lope and carrier representation of some field Q is 

Q{t-z) = yl(t;z)e*(*"'"^'^«*) -t- A*(<;z)e-'('^''^-"'"*),(l) 

In some of the following equations, I will shorten the ar- 
gument of the exponential in the carrier function with 
S = I {k^z — LUot). It is worth noting that we are not re- 
quired to use carrier functions with the usual exponential 



form [2^, e.g. in semiconductor physics, Bloch functions 
are routinely used as carriers to form a basis for electron 
(or hole) envelope functions. 

Note that it is approximations that restrict the valid- 
ity of envelope approaches, not the use of them in itself. 
This is contrary to the impression that might be gained 
from SVEA approaches, and even their generalizations 
17, 19]. The potential benefits of envelopes are not tied 
to restrictions on the bandwidth of the pulse being mod- 
eled. 



B. Tiie big advantage 

Replacing real fields E and H with an envelope-carrier 
description give us at least one clear advantage: it re- 
moves the dominant contribution to the underlying field 
oscillations. The resulting smoother envelope is there- 
fore easier, and much less computationally expensive to 
propagate. It is the rapidity of the fastest time-domain 
modulation of the field or envelope which constrains our 
time resolution, and the rapidity of the fastest spatial 
modulation which constrains the spatial resolution. Note 
that although we usually hope that our envelope will then 
have a relatively slowly varying form, the mere replace- 
ment of the EM fields by envelope-carrier combinations 
imposes of itself no approximations whatsoever. 

Two processes may act to twist an initially smooth 
envelope into something more problematic. First, linear 
dispersion can add chirp, which manifests itself as a non- 
linear phase ramp across the pulse. These are usually 
relatively smooth changes, and cause little problem. Sec- 
ond, there are nonlinear effects. Some, such as self phase 
modulation (SPM) can be relatively mild, others, such as 
coupling to backward propagating waves and harmonic 
generation can impose significant oscillations. 



C. Nonlinear polarization terms 

As mentioned above, nonlinear processes can generate 
oscillatory contributions to the envelope. We can see how 
this occurs by considering an instantaneous third order 
nonlinearity, which depend on and has the form 

E^{t;z) = [A(t; z)e- + A* {t; z)e--]^ (3) 
= Ait; zfe+^^ + 3A{t; zfA*{t; z)e+~ + c.((4) 

Here we see a useful side-effect of the envelope-carrier 
representation - that nonlinear terms can be separated 
into convenient components. In this example, the full 
X^^^ nonlinearity splits into a third harmonic generation 
(THG) term proportional to A^^ , and an SPM term pro- 
portional to \A\ A; along with complex conjugate coun- 
terparts (c.c). Clearly the THG term is non-resonant 
with the chosen carrier, and keeping such non-resonant 
nonlinear terms will impose significant oscillations onto 
our envelope as it propagates. Such oscillations break 
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approximations relying on a relatively smooth envelope, 
which is why in SVEA models they are discarded; how- 
ever there is no a priori requirement to do so. E.g., in 
the wideband Raman model of Kinsler and New [26, 2T\, 
the Stokes and anti-Stokes fields appeared as sidebands 
on the envelope spectrum. 



D. Multiple carriers, wideband envelopes 

We can generalize from a single envelope-carrier pair 
by using multiple envelopes with carriers at different car- 
rier frequencies. However, each added envelope greatly 
increases the number of individual polarization terms re- 
sulting from a nonlinearity, so as a rule it is best to 
use the minimum number possible. As an exercise, just 
calculate the number of terms in an eqn. ([3]) derived 
from a field defined as E = Aie^^^ + ^26+"^ + c.c, not 
E = ^ie+"i + c.c! 

Multiple carriers work best when there are multiple 
narrowband fields which resonantly interact, such as in 
an optical parametric amplifier [3, 0] (check) or for Ra- 
man processes [H, US] • In such cases we can ruthlessly 
discard nonlinear polarization terms which are not per- 
fectly in resonance with processes of our choosing. An- 
other use for multiple carriers is for including both for- 
ward and backward propagating fields. 

If we use multiple carriers, and also allow wideband 
envelopes, then it is possible for multiple envelopes to 
cover the same piece of spectrum. In a continuous math- 
ematical description this overlap will always happen, but 
in a discrete or numerical implementation it will depend 
on our parameters. 

This overlap is not necessarily a problem, as long 
as we are careful about assigning polarization terms to 
whichever envelope equation we choose - we must make 
sure not to add the same term twice, for example. This 
can also lead to a non-unique description, in the case 
where a polarization term could be equally well drive the 
evolution of one of two (or more) envelopes. Neverthe- 
less, such non-uniqueness does not break our model, it 
just allows us a choice which we might be able to use to 
our advantage. In a simulation, we could try to assist our 
numerics by managing the envelope spectra by generat- 
ing a total spectrum, and then reassigning components 
in the overlap region according to some smoothing pro- 
cedure. 

The use of wideband envelopes can raise some interest- 
ing issues. For example, if the bandwidth of an envelope 
is greater than its carrier frequency, then the envelope 
will extend into negative frequencies. This is not a prob- 
lem for our physical model, since we still have to recon- 
struct the field from the envelopes and carriers, and those 
negative frequency components are matched by comple- 
mentary positive frequency ones from the complex conju- 
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gate of the envelopC"'^ . Again, we might consider a spec- 
tral management scheme which swaps these unexpected 
components over, restoring the envelope to pure positive 
frequency content (and hence its conjugate to pure neg- 
ative frequency content). However, we then find that at 
zero- frequency we have introduce a hard cutoff in the en- 
velope spectra, and so induced unwanted oscillations in 
the time domain version of the envelope. 

However, while such spectral management might seem 
to offer advantages, in practice it makes little difference, 
and adds needless complication to simulation code. I 
would consider it only if some unexpected interaction 
was generating significant spectral content near the band 
edge of an envelope, at a position where it would be well 
within the spectral range of some other envelope; and 
even then it might be easier to simply increase the enve- 
lope bandwidth. 



E. Directionality 

A carrier imposes a direction of propagation, and most 
carrier-based models silently neglect even the possibility 
of backward propagating fields, even though there is a 
coupling between them. However, Casperson^] used 
both forward and backward carriers to construct an en- 
velope model with a separation of the forward and back- 
ward field components and interactions. The more recent 
paper of Sanborn et al. [291 used the same approach. 

Backward traveling components, if forced onto a for- 
ward traveling envelope, will appear as non-resonant 
terms. If identified correctly, these can then be discarded. 

See also section IVII for a discussion of the coupling 
between forward and backward waves which is induced 
by a nonlinearity. 

F. Moving frames 

In combination with a suitable moving frame, an en- 
velope representation can turn a pulse of rapidly oscil- 
lating fields moving at the speed of light into a smooth, 
nearly-stationary waveform - with commensurate gains 
in simulation speed. 

However, we need to guarantee that all contributions 
from backward traveling components are removed, oth- 
erwise the envelope will contain oscillatory components 
moving at approximately twice the frame speed. 

A typical moving frame is defined by for a frame speed 
V as t' — t — z/v. Thus the spatial derivatives in propa- 
gation equations are altered using 

dz = d^i - vdt (5) 



That is, they should be so matched. If they aren't, you've done 
something wrong. 
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G. Estimating the computational cost 

Consider a wideband pulse, with a bandwidth of the 
order of its centre frequency loq. In a full-field approach, 
this will have the fastest modulations of the field being 
of the order 2ujq. In comparison, an envelope approach 
results in the fastest modulations on the envelope be- 
ing of the order loq. Simplistically we might then hope 
that the envelope approach allows us to halve our time 
and space resolutions whilst still retaining numerical ac- 
curacy. For narrow band fields the advantage is much 
clearer - a bandwidth of wq/IOO might allow resolutions 
to be coarsened by a factor of 100. For fields of a wider 
bandwidth, we gain little advantage, unless we shift our 
carrier frequency luq to the centre of the spectrum, even 
if that is not co-incident with the dominant frequency 
component. A more comprehensive examination of the 
effects of numerical resolution has been given for the non- 
linear Schrodinger equation by Sinkin et al. 30]. 

Note that since the linear response (dispersion) of the 
medium can be done exactly in the frequency domain, re- 
gardless of step-size, it might seem more appropriate to 
focus more on the role of the nonlinear response when es- 
timating the necessary temporal and spatial resolutions. 
However the accuracy of a propagation method cannot be 
easily evaluated whilst ignoring the dispersive propaga- 
tion, because both dispersive and nonlinear effects occur 
simultaneously. Even if a split step method is used, they 
are interleaved, and their effects cannot be disentangled. 



H. Disadvantages 

The slight disadvantage of using envelopes is that real 
valued time dependent fields are replaced by complex 
valued envelopes. This doubles the amount of storage 
used during computations, and also requires the use of 
complex Fourier transforms rather than the faster real 
ones; although of course the spectra of the fields are com- 
plex in any case. In practice, the computational cost is 
small, although the complexity of the simulation code is 
increased. 



III. MAXWELL'S EQUATIONS 

When propagating fields in free space, we use the the 
source-free Maxwell's equations. To simplify the de- 
scription we transform their time-like behaviour into fre- 
quency space. This enables us to write the convolutions 
required to model the linear time-response of the medium 
(e.g. dispersion) as multiplications. However, since the 
form of the nonlinear response is not simplified by this 
process, a convolution in frequency space appears. In fre- 
quency space, time derivatives convert to factors of — loj, 
so the equations are 



dzEx{uj) = -iLofl{Lo')-k Hy{uj;z). (7) 
The denotes a convolution, 

g(T)*P(t) = [ Q{r)P{t-T)dT =J-^\q{u;)P{u;)\ .{8) 



A rather nice way to scale these equations is to define 
suitable e„ and /i„ corresponding to a suitably chosen 
refractive index, hence /i„ will usually be fiQ. This means 
c„ = l/e„M«, k = uj/cn, in = e{uj')/en /2„ = ji{uj')/nn- 
We then define e = ^ft^iE and h = ^/JI^H, which ensures 
e and h are of comparable sizes. This gives us the scaled 
Maxwell's equations 



dzhy{uj;z) 
dzCx {lo) 



-ikin{uj')-ke.j:{LL>;z), (9) 
-ikjln{uj') hy{uj; z). (10) 



It is worthwhile comparing this scaling with that from 
the directional fields approach in section IIVI with the 
correspondences y/e^ <-> and ^/Jmi f^r ■ 

For our purposes, there are two main ways to solve 
Maxwell's equations: either FDTDQ or Pseudo-Spectral 
Spatial Domain (PSSD)[i3]. In FDTD we propagate for- 
ward in time, holding the fields E{z), H{z) as a func- 
tion of space. However, in nonlinear optics, it is more 
convenient to use PSSD, where we propagate forward in 
space, holding the fields E{t),H{t) as a function of time. 
Under PSSD derivatives are calculated pseudospectrally 
[sH ]. However, its most important feature is that the 
entire time-history (and therefore frequency content) of 
the pulse is known at any point in space, so applying 
even arbitrary dispersion incurs no extra computational 
penalty. In contrast, FDTD (or other temporally propa- 
gated methods) must use convolutions or time-response 
models for dispersion. Although spatially propagated 
simulations (e.g. PSSD) make it difficult to incorpo- 
rate reflections properly, this is not a significant con- 
straint as most such simulations are only interested in 
uni-directional propagation anyway. 

For example, in a ID medium with linear dispersive 
properties defined by e^, /ir, containing a third order ;\;'^' 
nonlinearity defined by Cc, the equations are 

dzHy{uj;z) = -iujer{uj')Ex{u!;z) 

-iu{l,{uj')S[Ex{t-,zf] (c^)}*^,(c^;z)(ll) 
dzErc[uj;z) = -iujjlr{uj')Hy{uj\z), (12) 

where 'J[Q{t)]{uj) denotes the Fourier transform of some 
function Q{t). This model allows for the time-response 
of the nonlinearity, and is thus applicable to (weakly cou- 
pled) Raman systems as well. Note that the terms de- 
pendent on and are simple products. The linear 
dispersion combined with a time-dependent third order 
nonlinearity gives a permittivity function which would be 
written 



dzHy{uj; z) 



-iLje{Lj') -k Ex{(-o; z), 



(6) 



e(r,t) 
e(r, t)i.E(t) 



€r{T) + e,{T) k E{tf (13) 
e,(r) * E{t) + {e,(T) * E"" {t) ] E{i^A) 
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In the case of instantaneous nonlinearity, eE 
eeS'lE^]. 



irE 



A simple and efficient way to propagate these equa- 
tions is using staggered E and H fields, which allow us 
to use an Euler-like integration for each field, but achieves 
second-order accuracy 32]. However, while the E and H 
fields necessary for a forward propagating pulse are easy 
to determine for co- incident E and H, we need to use 
staggered initial conditions or else we get a significant 
backward propagating component. Even with correctly 
staggered initial conditions, we see a small spurious back- 
ward component, the size of which depends on the time 
step. This backward pulse is hard to get rid of com- 
pletely, but it can be filtered in the time domain when 
the two pulses have propagated apart far enough. An- 
other point to consider, particularly when generating the 
initial conditions for very short pulses, is the zero-force 
condition [33[. This can be easily satisfied by deriving 
the E and H fields for the pulse from a suitable vector 
potential, rather than simply assuming a form for the E 
field. 

When considering the solution of these these Maxwell's 
equations, it is useful to partly calculate the time deriva- 
tive of eqn. p4p . For dispersion and a time response x^^"^ 
this gives us three terms. 



dtie{T,t)*E{t)) 



dter{T)^E{t) 

+ {dt{e,{T)i.E\t)])E{t) 
+ {e,{r)^E^{t)]{d^E{t)llQ) 



Thus we see that to solve the equations, we will need to 
calculate the derivatives of three terms: the usual disper- 
sive term, the time response term, and the field. We will 
also need to retain the value of the time response term 
as well. Since the time response term contains a con- 
volution, it is best calculated in the frequency domain, 
which is particularly convenient when using pseudospec- 
tral derivatives. We will need two FFT's to transform 
E and E'^ into frequency space. There we construct the 
dispersion term and the time response term by simple 
multiplications, and set up arrays for the derivatives by 
multiplying by —lu. We then need four back transforms 
for a total of six in all: one more for the time response, 
and three for the time derivatives of the dispersion, time 
response, and field. For an instantaneous nonlinearity, 
we need only three FFT's: two forward transforms (for 
E and E'^), and one back transform for the combined 
derivative. In addition to these six (or three) FFT's 
needed to solve the dzH equation, the dzE equation re- 
quires another two, for a total of eight (or five). 



A. Envelopes 

Although it is not often done, we can represent 
Maxwell's equations using an envelope and carrier repre- 
sentation. We express the fields E and H using 

E.At;z) = A{t-z)e'''^''''-'^'''^ + A*{t-z)e'<^'''-'^'>%l) 
Hy(t-z) = F(t;z)e'('^«^-'^"*' +F*(i;z)e-*('=''^~'^''tl8) 

We insert these into the Maxwell's equations above, 
separate out the normal and complex conjugate (c.c.) 
parts, cancel the carrier exponentials present on both 
sides of the equations, and rearrange to leave only dz 
terms on the RHS, 

dzF{ijj;z) = -iuje{Lu')*A{LU]z)-ikoF{LU]z), (19) 
dzA{Lu;t) = ~iujil{uj')i^ F{uj;t) -ikoA{Lu;t). (20) 

Of course there is still much detail hidden in the per- 
mittivity e, since it contains the nonlinearity. Conse- 
quently, I do not apply this envelope definition to a 
general equation of motion because how e is expressed 
usually depends on the field and therefore on those en- 
velopes. Starting with eqn. (jl4p . and expanding e with 
terms for both (linear) dispersion and a time depen- 
dent x^^-* nonlinearity (ec) gives 

e(r,t) ★ A(t)e+" + c.c. 

= erir) i. {A{t)e+~ + A{t)*e~~} 
+ (ee(r) * {A(t)2e+2= 
+A{t)A{ty + A{t)*^e'^-}) 
X {A{t)e+~ + A{tyer~} (21) 

e{T,t)-k A{t)e+^ = e,.(r) * A(t)e+" 

+2[e,{T)*\A{t)f}A{t)e+^ 

+ {e,{r)i.A{tf}A{tre+~ 
+ {e,iT)i.Aitf}A{t)e+'^ (22) 
i{uj + UJq) A{uj) — er{u!' + UJq) -k A{uj) 

+2 {?,(cj' + cjo)J [|^(i)|'] iuj')^*A{Lj) 
+ {ic{oj' + u;o)3^[A{tf] {oj')}*A{u;r 

+ {iciuj' + 3c^o)y [A{tf] [uj')} ★ A{uj). 

(23) 



We can see in eqn. ([23)) that the first three of the terms 
(one dispersion and two SPM-like) are resonant with the 
chosen envelope, but the last (third harmonic generation) 
is not, and it modulates the envelope at 2a;o, (and subse- 
quently the propagation by ~ 2ko). Note the form of the 
second SPM-like term, which needs contributions from 
two A{t)'s and one A{t)* to have the correct frequency 
dependence, but the convolution is with the A{t)* and 
not an A{t) as might be expected. 
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Note that in the case of instantaneous x'"^\ the third 
RHS term reduces to Cc \A{t)\ A{t), giving 

€{t) i. A{t)e+~ = erir)* Ait)e+~ + See \A{t)f A{t)e+- 
+ecA{tfe+^^, (24) 
e{uj' + uJo) A{uj) — er(w + a;o)^(ijj) 



\A{t)fA[i) {uj) 



{A{tf\ 



(25) 



This expression can then be substituted directly into 
eqns. (fT9l^ 

In this formulation, we have made no "slowly varying" 
approximation like those in traditional approach es [ ij, \ A 
SSSIB], or in the variously corrected extensions jl7L[l3|. 
The price we pay is having two envelopes instead of one, 
since now the magnetic field is explicitly retained. Also, 
the model still contains backward propagating compo- 
nents; which, with the chosen carrier functions, will im- 
press oscillations at 2loq on the envelope, and oscillations 
of 2fco on the propagation, placing greater demands on 
our numerics. Unfortunately there is no way to filter 
these out at any point in the simulation, because their 
backward propagating nature can only be established by 
linking the time-like behaviour and space-like propaga- 
tion of both E and H fields. We cannot always rely on 
only time-like behaviour to filter them out, because, e.g., 
both backward propagating terms (at -|-fco and — o^o) and 
third harmonic generation (at -|-3fco and Swq) are equally 
detuned from the carrier (at -|-fco and wq); although we 
could do so if we were in a regime where third harmonic 
generation were negligible. 

At the start of this subsection, we hoped that dividing 
out the carrier oscillations would give us a slowly varying 
pulse envelope, which would then enable us to coarsen 
our numerical resolution, and speed simulations. This is 
true, up to a point - but remember the most likely reason 
we are using a Maxwell solver is that we want to model 
a wideband situation. It is the rapidity of the fastest 
time-domain modulation of the field or envelope which 
constrains our time resolution, and the rapidity of the 
fastest spatial modulation which constrains the spatial 
resolution. 

We can do better than these Maxwell equations ap- 
proaches without having to use second order wave equa- 
tions and their complicated approximations by using di- 
rectional Maxwell's equations, as described in the next 
section. 



B. Transverse effects 

There are two main transverse effect likely to be of 
interest in pulse propagation models: mode averaging, 
and diffraction or off-axis propagation. 

Mode averaging is easy to incorporate if you assume 
some known transverse profile for the mode: e.g. for 



an optical fibre or some other waveguide. The trans- 
verse derivatives vanish, and the material properties are 
evaluated as an integral over the transverse dimensions, 
weighted by the mode function. 

Diffraction and off-axis propagation they result from a 
coupling between the vector components of the E and El 
fields ~ including those along the propagation direction. 
Thus they are much harder to understand, as compared 
to a paraxial model based on (e.g.) the second order wave 
equation, although they can be simulated easily enough 
in a full 4D FDTD code. This is because they result 
from a coupling between the vector components of the 
E and El fields - including those along the propagation 
direction. 



IV. DIRECTIONAL MAXWELL'S EQUATIONS 

To my knowledge, the earliest rewriting of Maxwell's 
equations in a directional form was by Fleck [l^, who 
treated a dispersionless medium and plane polarized 
wave. However, the idea was not used beyond its brief 
appearance there. Fleck constructed his directional fields 
by combining the sum and difference of the E and H 
fields, weighted by the square roots of the permittivity 
e and permeability [i, respectively. The new combined 
fields represent the forward and backward traveling com- 
ponents of the total field, and we can derive first-order 
wave equations for these new fields. 

In the mid 1990's, the concept was rediscovered and 
used to evaluate the properties of grating structures by de 
Sterke, Sipe, and co-workers [s^fssj. but not applied to 
pulse propagation. The work considered materials with a 
spatially varying refractive index, but did not incorporate 
material dispersion or nonlinearity. 

This concept of using directional fields for pulse prop- 
aqation was not revisited until the work of Kolesik et al. 
pll IT^. After selecting a preferred direction, they then 
projected out the forward-like and backward-like parts 
of the propagating fields. This procedure resulted in first 
order wave equations for the propagation of the forward 
and backward field components. Subsequent work by 
Kinsler et al. [H, presented a directional rewrit- 
ing of Maxwell's equations using a generalized form of 
Fleck's construction. It is also worth noting the indepen- 
dent work of Mizuta et al. fl^, published at the same 
time. All of these methods use the same basic concept 
- use the right combination of E and H fields so as to 
create a pair of forward and backward-like fields. 

Here I follow the formulation of Kinsler et al. (l3l. 
which handles the electric and magnetic properties of the 
propagation medium on an equal footing, incorporates 
the dispersive properties of the medium in a very general 
way, and retains all the vectorial behaviour of the fields. 
The result is paired first-order equations for the plane- 
polarized directional fields (and a longitudinal com- 
ponent G°). Although complicated in the general case, 
these simplify greatly in the usual case(s) of transverse 
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and/or paraxial propagation regimes. The cost of us- 
ing these directional fields is that while we can efficiently 
remove backward propagating contributions, computing 
the nonlinear terms is more demanding. In contrast, the 
work of Kolesik et al. and Mizuta et al. is distinguished 
by a greater emphasis on the practical applications of 
directional fields. 

Because these new fields are directional, we can ef- 
ficiently separate out the forward-going part of the field, 
and neglect the backward. This is an important step, be- 
cause the standard Maxwell equations based approaches 
treated in the previous section could not easily remove 
the backward parts of the field, and these can cause in- 
convenience in numerical simulations. For example, the 
spurious backward component caused by imperfect initial 
conditions should no longer occur. 

The definitions of the fields, describing the trans- 
verse properties of a plane polarized EM field, in the 
frequency domain are 

G±(cj) ^ ar{Lo)E.Auj)±M'^)Hy{oj), (26) 

The ctr and Pr "reference" parameters are best chosen 
to closely match the medium, whilst ignoring nonlinear 
effects, so that ar{Lo)(3r{^) ~ ^/c{uj). That is, relevant 
(linear) dispersive properties of the medium are included 
in the reference parameters, i.e. that a.ri'^) = er('^)^^^- 
They have the definitions 

e ^ eriuo) + = a^(w) + Q:r(w) q;c(w), (27) 

jl = flr{uj) + flc{0j) = P^(uj)+(3r{uj) Pc{0j), (28) 

where the correction parameters Cc and p,c represent the 
discrepancy between the true values and the reference. 
These correction terms will usually just be the nonlin- 
earity. More generally, the smaller these correction terms 
are, the better the match, and the more likely it is that a 
description involving only G"*" will suffice. Note also that 
there are alternative ways of constructing directional G^- 
like fields PI. 

In the widely used moving frame defined hy v = 
l/af(3f, where dzQ = dz'Q — aff3fdtQ, using these Gj 
fields gives the (non-magnetic case) propagation equation 

El, 



dz>Gt = Tiw5,/3,(lT0 G± 



(29) 



where ^ — af(3f /urPr- Although this moving frame has 
no sensible limit as the frame speed tends to zero, the 
stationary frame case can be recovered by setting ^ = 
and replacing z' by z. G^ field simulations usually 
assume = 0, and treat only the forward traveling 
components of the EM field. 

Correctly writing down the form of nonlinear terms 
for eqn. (j29p requires some care, and consideration of 
the specific nonlinearity involved. Fortunately the task 



is simplified because it is simply a rewriting of the (elec- 
tric) nonlinear term from Maxwell's equations with the 
appropriate scaling factors relating ac to e, and Gj to 
E. 

Wave equations with a more familiar appearance can 
be obtained using 



E^iio) = G±(c.)/25.(c.). 



(30) 



These have the units of an electric field (i.e. V/m), but 
actually incorporate information about the magnetic field 
as well. If we take this step, we can transform back into 
forward propagating "electric fields" i?^, and get 

-dz'E^ = Ti^arPrilTd E^ 



E^ 



E' 



(31) 



An approximate forward-only wave equation can be 
found by setting E" = in eqn. ([?T|) . (or G^ = in 
eqn. (P^ ). For a time response x^^"* nonlinearity, this 
is 

-dz'E+{uj) ^ -lUjarPr {I - C) E+{lu) 

-iLO {/3.e~e(w).:r [E+it)^] iio)}*E+{u;) 

(32) 

Notice the similarity to eqn. pT|) . but that the field is 
propagated in a single first order equation, rather than 
two (i.e. both eqn. (fTTj) and p2|) ). The cost is that it only 
propagates forwards, but this is what we wanted. Fur- 
ther, the method can be implemented using fewer Fourier 
transforms than are required for a full Maxwell equa- 
tion solver The gain is that of not solving for dzE 
(eqn. (fT^ ). which requires a pair of FFT's if done pseu- 
dospectrally. Solving for pulse propagating in a medium 
with dispersion and a time dependent (or instantaneous) 
third order nonlinearity therefore requires only six (or 
three) FFT's, as compared to eight (or five) for solving 
Maxwell's equations. 

However, in practice the speed gain can be less clear 
cut. A PSSD solver moves forward one full step dz in two 
staggered steps, one integrating for the magnetic field, 
and integrating for the electric field; and only the mag- 
netic field integration needs to calculate the nonlinear- 
ity. This staggered scheme is second order accurate even 
though each stagger-step is only integrated using an Eu- 
ler method. We can achieve nearly the same level of ac- 
curacy for the directional fields by employing a leapfrog 
algorithm [3^ . If we wish to use more accurate (and so 
more complicated) numerical integration algorithms (e.g. 
a Runga-Kutta scheme), then we can only outperforms 
the staggered (or leapfrog) PSSD schemes if the propaga- 
tion step size is (greater than) twice that of the staggered 
PSSD. To complicate matters further, for reasons of nu- 
merical stability, we often need to tie the propagation 
step dz to the time grid step dt. This means that if there 
are bandwidth constraints limiting our dt, we may not 
have as much much freedom to adjust dz as we might 
like. 
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1. Special case: x'^' 

In the case of a x^^-* nonlinearity, two different field po- 
larizations are coupled together, and the equations given 
above tend to obscure the final form the nonlinear term 
will take. In this case, the time-domain displacement 
fields D in the two polarizations are 



€y -k Ey 



(2)£'2 



N(2)(,33) 
(34) 



If we assume that all of the linear response of the material 
(denoted above by ^XT^y) is absorbed into the reference 
parameters dr,f3r, we need only consider the nonlinear 
part. Note in particular that for the Dy field (i.e. J^y ') 
this does not depend on Ey, meaning that the forms of 
the wave equations given above (aimed largely at a x^^-* 
system) are not very useful. 

First, note that for a x*-^-* nonlinearity 



:^(3) 



(3)^2 



*E 



Ctr-Ctc * E, 



(35) 
(36) 
(37) 



and these give a nonlinear term for the wave equations 
of 



lUJOLcPr 



Gt + G- 



TUlPr.ardtc-k E, (38) 



By comparing these x^'^^ terms, we can see that in the 
X^^-* case, the nonlinear terms in the and wave 
equations (see eqn. (HH])) will be rewritten as follows 



X : ★ 



ILJOLcPr 

y ■ — n — * 



Gt ^ 


f G^ 




G+ - 


\-Gy 





2€oX^^^ExEy 



(2)£;2 



(39) 
(40) 



These can then be put in a form containing only G^, Gy 
if desired, but it is simplest to reconstruct the E^, Ey di- 
rectly before calculating the nonlinear terms. If a more 
extensive collection of the x^"^^ coefficients needs to be 
included, this procedure can be reproduced using the ap- 
propriate nonlinear field combinations. Further, if the 
time-response of the nonlinearity is also important, then 
we can include this by replacing x^'^^ExEy and x^'^^E'^ 
with appropriate convolutions: e.g. {x^"^^ * Ey)Ex and 

{x^-'UEx)Ex. 

For the £;±-like wave eqns. pip the nonlinear terms 
are 



Since we will want to apply the nonlinear effects in 
the time domain, we need to back-transform the terms 
in eqns. (|39|40p or eqns. (|41l42p . requiring a pair of 
Fourier transforms in addition to those required to get 
the time-domain fields. If using the G^ form, there is 
an additional transform, because we also need the time 
domain field(s) E{t) - with the E^ form E[t) can be 
found directly. 

Further simplifications can be made: e.g. in a semi- 
wideband limit around a central frequency wq, we can 
assume the frequency dependence of the a parameters in 
the nonlinear terms vanishes, so that the transform(s) to 
convert from G* to E is unecessary. In an SVEA-like 
narrowband limit all these transforms vanish because (in 
the nonlinear terms) the frequency dependence of the a 
parameters vanish and the the factor of uj simple becomes 
Wo- 



A. Envelopes 

Here I have intentionally simplified the definitions to 
best match what is most likely to be used in practice: a 
forward propagating G^ (or E^) only model. A more 
complete description of G^ envelopes, such as that in 
, would include the role of forward and backward trav- 
eling envelopes for both of G^ . 

We have seen that in the forward-only approximation, 
G+ and E^ follow identical equations of motion. The 
envelope and carrier representation of E^ is 

E+{t-z) = G(t;z)e'('=«^-'^''*) +G*(t;z)e-*('=''^-"'«t^43) 

I do not apply this envelope definition to the general 
equation of motion because how oic is expressed depends 
on the field and therefore on those envelopes. 

We now, for the case of a time response x^'^^ nonlinear- 
ity, substitute eqn. into eqn. ([5^ . we then (as usual) 
split the normal and c.c. parts, cancel exponentials, and 
rearrange leaving only the terms on the left, 

-zw/3,.ec(w-t-wo)-3' [2|G(t)|^j {uj).C{uj) 

- lUjPr£c{uJ + LUo).9 [C(tf] {lu).C*{lu) 

- iujPr£ciuJ + 3wo)-3' [G(i)^] (t^)-G(w), 

(44) 

The first line on the RHS will mostly cancel in the nar- 
rowband case, since /3,.ar — l/c(w), and kg — ll!q/c{ujo), 
thus with S ^ Lu — luq it becomes 



X : 


2 * 


'Et- 


^E- 




y ■ 


2 * 






25,. 



2eax^''^ExEy 
eoX^^El 



(41) 
(42) 



c(w) 



fco 



-I [k{u}) - fco] 
dk 



1 d^k 

2 aJ2 



(45) 
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where in the truncated expansion on the second line we 
can see the expected group velocity and group velocity 
dispersion terms. 

Note that eqn. is directly comparable to one de- 
rived from the NEE of Brabec and Krausz ^17|], but the 
only approximation I have made is to discard backward 
propagating fields. Since the NEE makes several addi- 
tional approximations, eqn. (j44p is more accurate and 
less approximate. Indeed, Brabec and Krausz were fortu- 
nate in that their chosen approximations produced a re- 
sult remarkably similar to that from the less restricted di- 
rectional fields approach. Note that Kolesik and Moloney 
[1^ also reduced their directional wave equation to a 
number of special cases, including that of Brabec and 
Krausz. 



B. Transverse effects 

As for Maxwell's equations, there are two main trans- 
verse effect likely to be of interest in pulse propagation 
models: mode averaging, and diffraction or off-axis prop- 
agation. 

Mode averaging is easy to incorporate if you assume 
some known transverse profile for the mode: e.g. for 
an optical fibre or some other waveguide. The trans- 
verse derivatives vanish, and the material properties are 
evaluated as an integral over the transverse dimensions, 
weighted by the mode function. This is just the same as 
for Maxwell's equations, although we now may be aver- 
aging slightly different quantities (e.g. ar rather than e). 
The work of Kolesik et al. [ll|, allows for transverse 
mode structure, that of Mizuta et al. [l^ for transverse 
averaging over a single mode. 

Diffraction and off-axis propagation are again much 
harder to understand, because (again) they result from 
a coupling between the vector components of the E and 
H fields - including those along the propagation direc- 
tion. However, second order wave equations derived from 
the first order directional fields equations exhibit a V'^ 
diffraction term which is the same as that seen in stan- 
dard second order wave equations (see e.g. eqn. (|47)) . 
in section |V|. This means that weakly transverse ef- 
fects can be accurately incorporated by using a split step 
scheme alternating between the wave equation and a 
diffraction term. Note that Kolesik et al. [l^l had wave 
equations incorporating diffraction (transverse) terms. 



V. SECOND ORDER WAVE EQUATIONS 

The standard second order wave equation applies to 
propagation in non-magnetic materials. If we consider 
the case of small transverse inhomogeneities of the po- 
larization, the three dimensional wave equation in typical 
notation (e.g. from [13, [3) is 



re Vj^ is the transverse Laplace operator, eL{t) — 
J^diJeL{u;)e^^\ e^ij) = I + 47rx(c^), and x(^) 



Here 

is the linear electric susceptibility. The electric field E 
propagates along the z direction. Both E and the nonlin- 
ear polarization P„; are polarized parallel to the x axis. 

Because of their starting point, methods based on this 
second order equation are slightly more restricted than 
those starting from Maxwell's equations. However, the 
differences in practice will likely be small, especially in 
the usual case of non-magnetic propagation media. 

Most uses of eqn. (HS)) . notably the slowly vary- 
ing envelope approximation (SVEA) relies on using an 
envelope-carrier description for the fields, then expands 
for weak dispersion, and resonant nonlinear perturba- 
tions about this carrier. This approach is discussed below 
in subsection IV Al 

Alternatively, we can attempt to factorise the equation 
into a product of two first order parts, as can be done for 
linear waves (see e.g. js^l)- Factorization is consider- 
ably more useful than the traditional approach, and is 
discussed below in subsection IV Dl 



A. Traditional approach 

Unlike the other approaches discussed in this paper, 
the traditional approach assumes the use of an envelope- 
carrier description of the field. 

Kinsler and New [I^, [^J] presented a comprehensive 
re-derivation of the envelope propagation equation based 
on the second order wave equation, which subsumes the 
SVEA and Brabec and Krausz's NEE [l3| as special 
cases. Since it is the most general, I use the Kinsler and 
New calculation, leaving some definitions to their paper 
rather than repeat them here. Noting that ^ and r are 
scaled space and time variables, that alpha and /? have 
different meanings from the rest of this paper, and that 
D' contains the dispersion terms, we have 

9^A(rl,e,r) 

^+^n')Air.,tr)+^jl^^Air,,tr) 
Po / (l + iadr) 



(47) 



where 
Tr 



A(fi,^,T).(48) 



{dl + Vl)E{r,t) 



dneL{T)^E{r,t)} 



An 
72" 



Eqn. (|47p is exact - it contains no more approxima- 
tions than the starting point eqn. except for the ex- 
pansion of e in powers of w. If we set Tr = 0, this gives us 
a generalized few cycle envelope (GFEA) equation, which 
contains the SVEA % Brabec and Krausz's NEE can 
be recovered from eqn. (j47p in the ID case where phase 
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and group velocities are the same (i.e. a = 1), likewise 
Porras's SEEA [l^ can be identified in the diffraction 
term. Of course we cannot just set the T/j term to zero 
without some justification, but this has already been ex- 
tensively discussed, not only in both but also the 
detailed analysis [2J]- 

Now consider the complicated few-cycle correction to 
the polarization term in eqn. (|47p . which contains par- 
tial derivatives (1 -I- lad-r) in the denominators. These 
will need to be evaluated by Fourier transforming into 
the conjugate frequency space (17). Further, the term 
is divided by another such term. Clearly these might, 
in wideband cases, result in denominators close to zero, 
causing the approximations to fail. This means they put 
a serious brake on the validity of any such approach, espe- 
cially if the bandwidth of the pulse approaches the carrier 
frequency. 

Note that the best first order expansion of the few- 
cycle corrections to the polarization term is more general 
than that given by Brabec and Krausz, and contains the 
group to phase velocity ratio cr, i.e. 



277r 



2m [1 + idr) 
nl (1 -t- ladr) 

Unfortunately for the venerable SVEA based on the 
second order wave equation, and even its most general 
variant presented here, the directional fields method dis- 
cussed in the previous section IIVI has made it utterly 
redundant; as, indeed, has the approach in the following 
subsection IV Dl There is no reason to use any form of 
the GFEA or SVEA when we can generate equations like 
eqns. (|31l32p by not only using fewer approximations, 
but much simpler ones than those taken by neglecting 
Tr. 



B. Time propagated direct solution 

It is of course possible to solve the second order wave 
equation by propagating it in time, either with or with- 
out the use of an envelope and carrier. This approach 
has been used with significant success by Scalora and co- 
workers (e.g. their early work [H, [H, Ho] ) . By propagat- 
ing in time reflections are handled correctly, an important 
feature when treating structured materials. Generally 
the solution is achieved retaining the second order spa- 
tial derivatives (both in z and transversely in x,y), but 
approximating the time derivatives to first-order. The 
approximation is made using an envelope with a well- 
chosen carrier frequency, and gives rise to the SVEAT, 
or slowly varying envelope approximation in time. 



C. Short pulse equation (SPE) 

The second order wave equation can be converted into 
the SPE by using a multiscale expansion ^J]- First, spe- 
cialize to a third-order nonlinearity (strength p) and only 



second order (ordinary) dispersion (strength d) and then 
rewrite the second order wave equation as 

1 



dtE{t; z) - —dtE{t; z) - d2E{t; z) ~ pdtEit; 



0(50) 



We introduce the scaled co-moving frame variables r = 
{t — z/c)/a so that dt = (l/cr)9r, and z„ = a"-z so that 
dz — —{l/cia)a'^dz^; hence after simplification eqn. ([50]) 
becomes 

--drdz,E{t-z)-d2E{t;z)~^dlE{t;zf = 0.(51) 

Ci <J^ 

Now, writing the field in multiscaled form as a power 
series in components Ei scaled by factors of cr, we have 

E{t;z) = (jEo{t,zi,Z2,...)+(j'^Ei{t,zi,Z2,...) + {52) 



and to leading order, we can write eqn. ()51l) down as the 
SPE 



-drdz.Ea - d2Ea - pdiE^ = 0. 

^1 



(53) 



This equation has spawned a literature all of its own, 
because it (like the ordinary nonlinear Schrodinger equa- 
tion) provides a rich variety of mathematical solutions. 



D. Factorization approach 

An alternative to the traditional style of derivation dis- 
cussed above, we can instead factorise the second order 
wave equation in a way similar to that done for linear 
waves (see e.g. [13]). This was initially suggested by 
Shen [2], followed by Blow and Wood i2u\, and more re- 
cently revisited by Ferrando et al. and Genty et al. 

First we reduce eqn. (j46|) to the ID strictly paraxial 
limit; then transform into frequency space. Here I re-use 
the symbol (3 as the propagation wave vector to match the 
notation of Genty et al. [23], so that /3(w) = LUy/er{uj)fio. 
The wave equation therefore is 

5f i?(t; z) - \d^E{t; z) - ^lo^!Pit■, z) 



0,(54) 

0,(55) 
0,(56) 



(3) 



V^E{uj; z) + P{ujfE{uj; z) + fiQUj'^P{uj; z) 
dlE{uj; z) + {lj)E{uj; z) + (w)DSf ★ E{uj; z) 
where for a third order nonlinearity, with ec = eoX 

li - fioeoX^''^uj^E{r,t;zf/p{u;f 
= fioeci^^Eir,t;zf/Piu;f 

^ JL^E{r,t-z)\ 

I now briefiy consider three factorization approaches, 
from the simple method of Blow and Wood [20] , an im- 
proved version, and finally the most rigorous approach. 
Although these traditionally involve an envelope-carrier 
decomposition introduced early in that calculation (see 
Blow and Wood), the step is in fact unnecessary and I 
omit it. 



(57) 
(58) 

(59) 
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1. Simple factorization 



Linear factorization 



Factorization approaches are simple in two situations: 
a dispersionless medium with an instantaneous nonlin- 
earity, and a dispersive medium with no nonlinearity. In 
the dispersionless nonlinearity case, we can factorise in 
the time domain. In the linear dispersive case, we can 
factorise in the frequency domain. In the dispersive non- 
linear case, it is (usually) not possible to analytically fac- 
torise the second order wave equation. 

Regardless of mathematical difficulties. Blow and 
Wood [l^l used a simple factorization procedure which 
ignored the details of nonlinearity and dispersion. Re- 
membering that f3 — f3{uj), and without their envelope- 
carrier decomposition, they had 



E = 0. (60) 



They then separated out the forward propagating term. 
The envelope equivalent of this was then expanded us- 
ing a "weak nonlinearity" assumption with a binomial 
expansion, keeping only the first order corrections. 



2. Improved factorization 

The approach of Blow and Wood ignores the mathe- 
matical difficulties due to the use of the square root in 
combination with the frequency-domain convolutions be- 
tween the nonlinear term DnT and the field spectrum E; 

Fortunately, we can instead "complete the square" 
(e.g. l + iV~ 1 + A^ + = (1 + N/2f), enabling 

us to preserve the convolutions correctly. This requires 
us to make a weak nonlinearity approximation, but it 
is one nearly identical to that used when expanding the 
square root in the Blow and Wood calculation. So, with 
a the weak nonlinearity constraint 



^Sii.E<S:l, 



we get 




E 



(61) 



062) 



By assuming the forward-like and backward-like terms 
in square brackets factorise. 



E 



(63) 



d^E = ±il3E ± ip'^E. (64) 

While this equation can give excellent results, it is re- 
stricted to weak nonlinearity: as we see below, it lacks the 
nonlinear coupling term between the forward and back- 
ward propagating fields. 



Ferrando et al. [2l| treat this approach in detail, sep- 
arating this second order equation into two first order 
equations using a Greens function approach. More re- 
cently, Genty et al. [l^l have used it to generate a non- 
linear envelope equation, and specialised it to the case of 
nonlinear waveguides. First we move the nonlinear term 
from eqn. (|56p to the RHS and Fourier transform z into 
fc-space, so that the dz becomes —ik, 



[-k^ + l3^]E 
E 



-(3^Ji-kE 



yi-kE 



{k + /3) (fc - /?) 



E++E^ = 



2/3 



1 



k + (3 k-(3 











E+ 


+ E^ 



(65) 
(66) 

(67) 
(68) 



where I've taken i? to be a sum of both forward and back- 
ward propagating parts E = E^ + E- . I now split eqn. 



into a sum of two parts, where each half represents 
the propagation of the forward field E+ or the backward 
field E_, and rearrange. 



E. 



E+ + E^ 
+ E^ 



Now I transform back from /c-space into z, so that 



E+ + E^ 



E++E^ 



(69) 
(70) 

(71) 
.(72) 



If we compare this result (i.e. eqn. (ff^ ) with the 
comparable equations for the directional fields G^, in 
particular with the electric field form given in eqn. (j3ip : 
we see that they are essentially identical: since LuarPr — 

Uj/Cr ^ (3. 

A similar procedure can be applied to eqn. (|62p if 
desired. Eqn. (I72p is almost the same as the (more ap- 
proximate) eqn. (l63|) . Since the RHS nonlinear term is 
a function of {E+ + ES), it provides a route for coupling 
between the forward and backward waves; its form can 
be obtained from the nonlinear part of (e.g.) eqn (|lip . 
A specific example for the case of a time-response x*-^-* 
nonlinearity has been given in |22l | , but in my notation it 
is identical to that for the rescaled directional fields 
(i.e. £;±), i.e. eqn. ([5^ . 



4- Special case: x'^' 

In the case of a x*^^'' nonlinearity, two different field 
polarizations are coupled together, and the equations 
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given above tend to obscure the final form the nonlin- 
ear term will take. First, note that the factorisation pro- 
cess changes the nonlinear term from (j'^uP'yi * E into 
iPuJ^'H -k E /2. This means that the term itself is multi- 
pled by a factor of just «/2/3, and this transforming factor 
is what we need to use in the general case. 

For a x*-^^ nonlinearity, the time-domain displacement 
fields D in the two polarizations are 

D,, = e.,*E^ + 2eQX^^^E^Ey, (73) 
Dy = ey*Ey + eoX^^^El (74) 

Note in particular that the nonlinear part of the Dy field 
does not depend on Ey, making the wave eqn. ((72|) 
(aimed largely at a x'^' system) inappropriate. 

In any case, the x nonlinear term is just 2eox'"^^ E^Ey, 
and the y term eoX^^''^^ so that in the pair of frequency 
domain wave equations (cf eqn. (j55p ). the nonlinear 
terms are 



y 



and in the factorised equations these become 

, ,2 



X : 



y ■■ 



/3(c.) 



X^^'^E^Ey 



(75) 
(76) 

(77) 
(78) 



Since we will want to apply the nonlinear effects in the 
time domain, we need to back-transform these nonlinear 
terms: 

, ,2 



X : 



y ■ 



Pico) 



UJ 



2(3{lo) 



X 



^^^E^Ey 



(2) 



Ei 



(79) 
(80) 



So we see that a true wideband approach to the non- 
linearity requires a pair of Fourier transforms. In a semi- 
wideband limit around a central frequency loq we can 
probably assume the factor Lu'^/f3(uj) becomes cuj / n(ujQ) . 
In an SVEA-like narrowband limit it would become 
Wq//3(wo) — cujQ / 'n{ujQ) , and the need for Fourier trans- 
forms vanishes. 

If a more extensive collection of the x*^^"* coefficients 
needs to be included, this procedure can be reproduced 
using the appropriate nonlinear field conbinations. If the 
time-response of the nonlinearity is also important, then 
we can include this by replacing x^'^^E^Ey and x'^^^E'^ 
with appropriate convolutions: i.e. (x''^^' * Ey)Ex and 

{x^-'Ue,)e,. 



tion is linear in the derivatives, when split into and 
A*_^ parts it looks very similar, being 



A, 



(81) 



For the case of a time-response x'^'' nonlinearity, the 
equation will be identical to that for the envelope version 
of the directional fields, i.e. eqn. 



6. Factorized fields 

An important feature of this approach is that we see 
that any contribution (whether linear or not) that is in- 
cluded in the source term will couple the forward and 
backward fields together. Consider two differing factori- 
sations of the same systems; e.g. one with the loss in- 
cluded in /3, and one with it in the source term. The one 
with the extra source contribution will see a correspond- 
ing extra forward backward coupling term, apparently 
conflicting with the fact that the two factorisations are 
of the same system. The resolution of this conundrum is 
simply that the forward and backward fields of the first 
factorisations {Ei±) are not the same as those for the sec- 
ond (i?2±); the meaning of "forward field" (or "backward 
field") differs between the two implementations. This is 
perhaps clearer in the formulation (see section ITV| . 
where the different factorisations would correspond to 
different choices of the reference parameters , /3r • If no 
further approximations have been made, when the real 
electric and magnetic fields are reconstructed from any 
factorised Ei±, the answers should be in agreement. 



E. Transverse effects 

In common with most pulse propagation, we can 
restrict ourselves to paraxial beams and incorporated 
transverse effects by using a split step scheme alternat- 
ing between the wave equation and the diffraction 
term. This is equally applicable to either the traditional 
or factorization approaches. However, in the factoriza- 
tion approach we can treat the V^^i? diffraction term as 
a "source" term, and, like the nonlinearity, move it to 
the RHS before factorising. Thus eqn. (|72p could be 
rewritten to include diffraction as 



id ~ 

d,E± = ±t(3E± ± yN* 



E. 



± 



2/3 



E. 



E^ 



(82) 



VI. FORWARD-BACKWARD COUPLING 



5. Factorization and envelopes 



Taking only the forward part of eqn ([72|) . we replace 

i?_i_(ti;) — A^{lu + loq) + A*,{uj — luq). Since this the equa- 



We can see in eans. ([29l [3T1 [72t that we simplify into a 
forward-only picture by dropping the part of the nonlin- 
ear polarization term due to the backward field. In situ- 
ations where there is no pre-existing backward field, and 
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where there are no interfaces to cause reflection, this is 
an excellent approximation that holds true in the regime 
of weak nonlinearity. It is only an approximation, be- 
cause the nonlinear polarization drives both the forward 
and backward fields, so in strongly nonlinear systems, 
a backward wave can be generated directly by the for- 
ward wave. The important "weak nonlinearity" criteria 
for perturbative nonlinearities to guarantee the validity 
of a forward-only model is 

^ J2 mx^"'^E"'-^ < 1. (83) 

'^0 m>l 

On the subject of reflections from interfaces, it is worth 
noting that "nonlinear" reflections can occur even if the 
linear dispersion on both sides is identical - as long as the 
nonlinearity changes, as in e.g. periodic poling, where its 
sign changes. 

Note that since nonlinearities are in practice very weak 
(e.g. ~ 0.06 at the damage threshold of fused 

silica), uni-directional propagation models perform very 
well, and the role of nonlinear reflections is generally neg- 
ligible. 

VII. CONCLUSIONS 

I have described three forms for the spatial propaga- 
tion of optical fields: Maxwell's equations, directional 
fields, and second order wave equations. These forms 
have been describe in both standard and envelope-carrier 
pictures. While solving Maxwell's equations remains the 
"gold standard" and most exact procedure, it is compu- 
tationally demanding, and it can be difficult to set up 
initial conditions. These difficulties are avoided by using 
a directional fields approach, where we can propagate 
more efficiently in the usual forward-only cases. Further, 
envelope theories based on forward-only directional fields 
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give equations of motion similar in form to the traditional 
SVEA ones based on the second order wave equation, but 
without requiring complicated approximations. 

When comparing the various approaches taken to di- 
rectional fields, a number of important points stand out. 

1. The first successful attempt at deriving useful di- 
rectional versions of Maxwell's equations was by 
Kolesik et al. [ll[i2. 

2. The most flexible and complete formulation is the 
directional fields of Kinsler et al. [l^, relying 
only on simple combinations of Maxwell's equations 
to achieve a directional form. It is applicable to 
propagation media with any frequency-dependent 
electric or magnetic properties, and variant forms 
(l^ can be used if required. 

3. The factorization style approach [13, HH, H^l gives 
propagation equations for the electric field that can 
be simply expressed and solved for without the con- 
struction of the somewhat opaque fields used by 
Kinsler et al. 

It is encouraging that these three approaches discussed 
in this paper (Maxwell's equations, directional G^ fields, 
and factorized second order wave equations) all give es- 
sentially identical results in the case of uni-directional 
propagation in non-magnetic media. 
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